Noise filtering method for seismic data

ABSTRACT

A method for filtering noise from discrete noisy seismic signals in a time interval ([1, . . . , T]) is provided. The method comprising the steps of using the noisy seismic signals for determining at least one reference channel (x(t)) for the time interval as an estimate of the noise; determining coefficients for M temporally local filters (w(i, t)), the filters forming a filter bank, and M being a number equal to or larger than two, by minimizing a cost function (J(t)) representing a measure of the error of the output of the filter bank, and applying the filter bank to at least one estimate (x(t)) to determine M filtered estimates of the noise. The filtered estimates are multipled with temporal windows (h(i, t)).

The present invention relates to a method of filtering noise from datasignals, and in particular to the filtering of noise so as reduce theeffect of coherent noise on data signals acquired during seismicinvestigations.

BACKGROUND OF THE INVENTION

In land seismic surveys, a source induces seismic waves at or near thesurface of the earth. These waves propagate through the earth andreflections from different layers within the earth can be detected bysensors, or geophones, at the earth's surface. The seismic sourcevibrations applied to the earth's surface also generate a so-calledsurface wave or ground roll which propagates through the shallow layersof the earth. At the geophones, the time of incidence of the lowfrequency, low speed ground roll typically coincides with the incidenceof reflections from the deep layers of interest in the seismic survey.The simultaneous presence of the ground roll with the reflected signalsmakes it difficult to make full use of the seismic data as the groundroll often masks the reflected waves.

Several methods are known for attenuating ground roll interference andthus reducing its effect on the seismic signal of interest. Typically,geophones are not used individually, but rather are connected insub-arrays, or groups, which are hard-wired or summed together. This isa form of data-independent beamforming. Attempts have also been made toapply adaptive signal processing for the suppression of ground roll inseismic surveys.

U.S. Pat. No. 4,556,962 attempts to attenuate the ground roll from asurface seismic source by placing a sensor close to the source to detectthe interfering noise. The interfering noise is scaled, delayed andsummed with signals from a more distant geophone array and thencross-correlated with the original vibrational source. This Patent alsosuggests that an adaptive filter be used so as to modify the delayedsignal to correspond more closely to that detected by the more distantgeophone array. However, the ground roll measured close to the sourcemay be substantially different from that received by the geophone array,and the adaptive filter may not be able to deal with this.

In U.S. Pat. No. 4,890,264 a method for suppressing non-uniformlydistributed noise generated by surface wave propagation, wind, andmachinery is described. Horizontal geophones for detecting surface wavesare used with conventional vertically orientated geophones for detectingseismic energy. The outputs of the surface wave detectors are used inconjunction with an adaptive filter to cancel the effects of the surfacewave interference. This method for the suppression of ground roll isinherently a multicomponent method. Some seismic wave energy also getsdetected by the horizontally sensitive geophones, and this may causesignal cancellation.

In UK Patent Application GB-A-2273358 the use of linearly constrainedadaptive beamforming and adaptive interference cancelling beamformingfor ground roll suppression was proposed. This method filters signalsmeasured by an array of geophones and sums them in such a way as topreserve signals incident from a preferred direction while suppressinginterference incident from other directions. The filtering is performedusing a continuously adaptive method with the moveout differentialbetween the seismic reflections and the ground roll being used to formprimary and reference channels. The suggested application is in drillingwhen using a drill as a seismic source, where the ground roll iseffectively stationary due to the slow travel of the drill bit and eachsource receiver position produces a lot of data. This ensures that thestochastic gradient type of algorithms used in the adaptive filters ofthis method are able to converge. However, in surface seismicexperiments the ground roll present is often non-stationary andinhomogeneous and the stochastic gradient type of algorithms may be tooslow to converge within the signal envelope.

U.S. Pat. No. 5,237,538 proposes a method for removing coherent noisefrom seismic data. This method firstly identifies the moveoutcharacteristics of the noise, defines and extracts a space-time gatecontaining the noise, and removes the moveout to flatten the noisetrain. Amplitudes and time variations are then removed from the gate.The coherent noise is estimated using a beamsteer operator (conventionalstacking in this case) or by f-x filtering in the Fourier tranformdomain. The filtered noise estimate is subtracted from the data tracecontaining the signal-plus-noise using a short three to five pointsingle filter. Inverse amplitude scalars are applied to undo the effectof earlier amplitude equalisation. The signal is then moveout restoredinto the original seismic record. This method has some particularshortcomings for application for ground roll attenuation. First of all,especially for shorter arrays, the signal always leaks into the groundroll estimate. In fact, there is always a component of the signalpresent at the reference channel which is colocated in time with thesignal in the primary channel. On the other hand, when the arrays areallowed to be longer, the dispersion present in the ground roll make itdifficult to achieve effective beamsteering.

SUMMARY OF THE INVENTION

In accordance with the present invention, a method for filtering noisefrom discrete noisy seismic signals in a time interval ([1, . . . , T]),said method comprising the steps of

using said noisy seismic signals for determining at least one referencechannel (x(t)) for said time interval as an estimate of said noise;

determining coefficients for M temporally local filters (w(i,t)), saidfilters forming a filter bank, and M being a number equal to or largerthan two, by minimizing a cost function (J(t)) representing a measure ofthe error of the output of said filter bank; and

applying said filter bank to said at least one estimate (x(t)) todetermine M filtered estimates of said noise.

Preferably the cost function (J(t)) is temporally global, representing ameasure of the error of the total output of said filter bank in the timeinterval([1, . . . , T]).

Further the cost function (J(t)) may be minimized using theapproximation that the sum, weighted by window functions, of the outputof adjacent filters of the M filters (w(i,t)) is equal when applied tothe same signal in time regions where said window functions overlap.

Preferably the method includes the step of multiplying M filteredestimates with temporal window functions ((h(i,t)).

The temporal window functions may be characterized by the requiementthat only adjoining windows overlap.

The application of the temporal window functions, and hence theresulting temporal windows, to the combined components ensures that thefiltering means is local in time and allows the method adaptively toremove noise from the seismic data in accordance with a globaloptimisation criterion, e.g. to solve the optimisation of the filteredsignal by minimising the mean square value of the filtered signal overtime.

The invention is applicable for two-dimensional (2D) andthree-dimensional (3D) seismic surveys, and can be used in land seismic,marine seismic including sea bottom seismic, and transitional zoneseismic.

The method can be performed on stored data or on raw seismic data as itis acquired. Thus raw seismic data may be filtered according to themethod at the data acquisition site. This ensures that a "cleaned"signal is available from the data acquisition site and may be downloadeddirectly from the site in this form. This reduces the amount of datathat must be sent for analysis off-site and reduces the costs andstorage problems associated with accumulating sufficient quantities ofnoisy data for analysis off-site. The method can be applied tosingle-sensor recordings, i.e. to recordings prior to any group formingwhich combines the signals of two or more seismic sensors.

The noise or reference signal may be pre-processed before being passedto the adaptive filtering means by dividing the signal into differentfrequency bands, for example by using a quadrature mirror filter. Thisallows a reduction in the number of data points to be processed and alsoallows a reduction in the number of coefficients in the adaptivefiltering means.

The data selection temporal window functions are preferably determinedby two requirements, wherein the first requirement is that the sum overall windows at any given time equals unity and the second requirement isthat only adjoining windows overlap. These requirements ensure that theglobal optimisation of the filtered signal can be solved by use of anapproximation in which for the sum over all time and all filters and allneighbouring filters, the error function of a neighbouring filter isreplaced with the error function associated with the filter itself.

The application of the data selection temporal windows decouples theequation required to solve the optimisation of the filtered signal. Theresulting decoupled equation may for example be solved by the method ofprincipal components or alternatively a damped least squares approach.

Where the decoupled equation is solved using the method of principalcomponents, the number of principal components may be adjusted to varythe degree of filtering and achieve the desired accuracy of filtering.The smaller the signal to noise ratio, the greater the number ofprincipal components that are used to achieve filtering.

Preferably the adaptive filtering is achieved by use of a filteringmeans, or filtering bank, comprising a plurality of local filters. Eachlocal filter may be multichannel filter. In a further preferredembodiment of the invention, the another signal is used to partiallyconfigure the adaptive filtering.

In land based surveys, the seismic data signals for use in methods inaccordance with the invention may be acquired from at least two seismicsensing means, typically geophones, the sensing means being arranged ina net, where each sensing means is independent of the other, and eachsensing means sends data signals to a processing means on site to befiltered in accordance with the invention. A net is defined as an arealarray of sensing means, or geophones, where each sensing means isindependent of all other sensing means in the net. This differs fromprior art arrangements where each geophone is arranged in aninterconnected array and the signals received by the geophones areaveraged over the array in an attempt to reduce noise effects such asground roll prior to data processing off-site.

In a further embodiment of the invention, the filtered output signal maybe reprocessed in an iterative manner to a further filter the noise, thefilter output signal typically being fed back to at least one of thereference signal and filtering means.

The reference signal may be generated by a number of techniques, ofwhich one example is by moveout differentiation of the data signals.

In another technique the reference signal may be obtained by medianstacking. This suppresses seismic signals that do not have the samemoveout as the noise to be filtered. Thus the contamination of thereference signal with the seismic signal of interest is reduced.

In a combination of these two techniques, median stacking may befollowed by moveout differentiation of the stacked signals.

Where the seismic data is obtained from multicomponent sensing means,the reference signal may be obtained by polarisation filtering of eachcomponent sensed by the sensing means. This method may be used toenhance the reference signals from adjacent sensing means which is ofparticular advantage for three dimensional exploration, i.e. threecomponent, and may also be combined with moveout differentiation.

These and other features of the invention, preferred embodiments andvariants thereof, possible applications and advantages will becomeappreciated and understood by those skilled in the art from the detaileddescription and drawings following hereinbelow.

DRAWINGS

FIG. 1 shows a processing block diagram of a prior art method foradaptive cancelling using a beamforming method;

FIG. 2A shows a method of adaptive interference cancelling in accordancewith the invention;

FIG. 2B shows a variant of the method of FIG. 2;

FIG. 3 shows sampling of the data signals to support the filtercoefficients used in the adaptive interference cancelling;

FIG. 4 shows a data-adaptive multichannel filter bank according to theinvention applying moveout differentiation to the reference channel;

FIG. 5 shows an alternative sampling of the data signals to support thefilter coefficients used in the adaptive interference cancelling;

FIG. 6 shows generation of the reference channel based on moveoutdifferentiation and median stack in the method according to theinvention;

FIG. 7 shows generation of the reference channel based on median stackand moveout differentiation;

FIG. 8 shows generation of the reference channel based on polarisationdifferentiation; and

FIG. 9 shows generation of the reference channel based on polarisationand moveout differentiation.

EXAMPLE(S)

Referring now to the drawings, FIG. 1 shows a linearly constrainedadaptive multichannel beamforming array as described for example inGB-A-2273358. The system inputs 10, 12 are derived from two sets ofsensors, one input 10 being from a primary source or sensor andsuccessive inputs 12', 12", . . . , 12^(N) from reference sources orauxiliary sensors. The input 10 comprises an information bearing signalcorrupted by noise.

The reference inputs 12', 12", . . . 12^(N) are processed by an adaptivemultichannel filter comprising local adaptive filters 14', 14", . . . ,14^(N) and the processed components combined at a summing means 16 toproduce an output signal y(t) 18.

The filtering is performed using a time-continuously adaptive method,i.e. the least-mean-square (LMS) method to optimize the performance ofthe filter. However, stochastic gradient type algorithms such as leastmean square may be too slow in converging to produce an estimate of thesignal without noise when incoherent noise is non-stationery andinhomogeneous, such as those encountered in surface seismics. Groundroll noise is in general: dispersive, and as such often displaysnon-stationary behaviour.

The method of adaptive interference cancelling according to theinvention estimates the phase and amplitude perturbation effectsalongside the seismic signal propagation effects while the adaptivefilter coefficients are computed. This makes the adaptive interferencecancelling approach robust to phase and amplitude perturbations presentin the data.

FIG. 2A shows a system for performing a method of filtering noise from aseismic signal in accordance with the invention. The system comprises aplurality of geophones G_(n), n=1, . . . ,N 24', 24", . . . 24^(N),which record signal of the general form {g_(nj) (t): n=1, . . . , N;j=1,2,3}, where j refers to a multicomponent recording, e.g a 3Cgeophone recording. To avoid confusion with FIG. 1, it should be notedthat the letter N now is used to denote the number of sensors, whereasthe letter K (see below) refers to the maximum number of referencechannels.

The geophones 24 are separate sensors placed on the earth's surface in anet to detect seismic waves produced by reflection at layers within theearth in response to a source wave. The geophones are not linked to eachother and act to provide a number of inputs to the primary and referencechannel generator 26. Equivalent systems using pressure sensors, e.g.hydrophones, instead of geophones may be used.

A signal generator 26 receiving said recorded signals designates primarysignal channel d(t) 28 which carries the information-bearing signal s(t)corrupted by additive interference n₀ (t) according to

    d(t)=s(t)+n.sub.0 (t)                                      [1]

and interference reference signal channels x_(k) (t) 30:

    x.sub.k (t)=F{g.sub.nj (t)}                                [2]

with k=1, . . . , K denoting the number of reference channels. Theoperator F can be any operation acting on the recorded signals so as togenerate one or more reference channels, which can be regarded as firstestimates of the additive interference. Hence, it is assumed for thefollowing that the correlation of the signal s(t) with the interferencein the primary channel, n_(o) (t), and with the reference channels x_(k)(t) is negligible compared to the correlation of n_(o) (t) with x_(k)(t). Examples for the operator F can be chosen from a broad variety ofmethods known as such. Some particularly advantageous methods forgenerating the reference signal 30 are discussed below.

The reference channels x_(k) (t) are processed by a number of localadaptive filters 32', 32", . . . 32^(k) arranged as an adaptive filterbank 34, summation means 36', 36", . . . 36^(M) , multiplication means38', 38", . . . 38^(k) and combining means 40 to produce the outputsignal ##EQU1## where w_(ik)τ (t) are the adjustable coefficients of theadaptive filters, h_(i) (t) are the windows applied at the output end, Mis the number of local multichannel adaptive filters (or the number ofoutput windows), and L is the number of coefficients per channel. Eachadaptive filter may be a multichannel filter of the type shown in FIG.1.

Equation [3] can be rewritten as a (windowed) sum over a scalar productusing a tap-input vector x(t) at time t defined as:

    x(t).tbd.[x.sub.1 (t), . . . , x.sub.1 (t-L+1),

    x.sub.2 (t), . . . , x.sub.2 (t-L+1),                      [4]

    x.sub.K (t), . . . , x.sub.K (t-L+1)].sup.T

and a tap-weight vector w_(i) defined as

    w.sub.i (t).tbd.[w.sub.i10 (t), . . . , w.sub.i1(L-1) (t-L+1),

    w.sub.i20 (t), . . . ,w.sub.i2(L-1) (t-L+1),               [5]

    w.sub.iK0 (t), . . . , w.sub.iK(L-1) (t-L+1)].sup.T

The filter bank output is combined with the primary input to produce afiltered output signal e(t) 44, where

    e(t)=d(t)-y(t).                                            [6]

In FIG. 2A the primary input is summed at the output of each filter 32.An equivalent method for determining e(t) is illustrated by FIG. 2B,where the window function h_(i) (t) is multiplied with the output ofeach filter, and the primary channel is added after all outputs arecombined.

The signal e(t) is also known as the error signal and a cost function Jis defined which is equal to the sum of e² (t) over all the sampleinputs i.e. ##EQU2## where T is the total number of samples. The costfunction J corresponds to the total error energy and by minimising thecost function, and thus minimising the global criterion, e(t) becomesthe least square estimate of s(t), i.e. the pure signal with no noise.

One way of minimising the global criterion is by using the windows h_(i)(t) to ensure that the filters are local in time, i.e. the arrays have afinite time span, and can deal with non-stationary noise. The windowsmay be triangular with a length of 100 samples.

The windows 42 are constrained by ##EQU3## for t=1, 2, . . . , T, whereT is the total number of samples, and

    h.sub.i (t)h.sub.j (t)=0                                   [9]

for j≠i-1,i,i+1. The first constraint ensures that the filter bank 34 isequivalent to the single filter case if all the local filters 32 areidentical. The second constraint ensures that the windows have compactsupport. Use of these two constraints enables the global optimisation ofthe error signal e(t) to be solved.

The optimisation problem is ##EQU4##

The total error energy can be expanded as ##EQU5## where f_(i) (t) isthe local error function associated with local filter w_(i) :

    f.sub.i (t)=d(t)-w.sub.i.sup.T (t)x(t).                    [12]

It can be seen from equation [11] that the optimisation problem leads toa large number of coupled equations for the filters. Using the secondcondition, eq. [9], imposed on the windows and the followingapproximation [13] for the second term in equation [11] the equationscan be decoupled. ##EQU6##

This is a very mild approximation as firstly, the approximated secondterm in [11] is dominated by the first term and secondly, theapproximation is milder than the approximation w_(i) =w_(i+1) or w_(i)=w_(i-1). The approximation of equation [13] requires that neighbouringfilters produce similar results when applied to the same input data intime regions where adjacent windows overlap, instead of requiring thatneighbouring filters are similar on a point by point basis. Thus, theapproximation is similar to requiring that the integral of two functionsare close, rather than the functions themselves.

Whilst optimisation leads to a large number of coupled equations for thefilters, the approximation of equation [13] decouples the optimisationso that ##EQU7## with i=1,2, . . . M.

The matrix normal equations for the components of the filter bank can beobtained, for example by setting the gradient to zero, or by using theprinciple of orthogonality, to give: ##EQU8##

This can be solved by various methods, two of which are the method ofprincipal components and the damped least squares approach.

The right hand side of quation [15] can be interpreted as a product of acorrelation matrix Φ and the filter coefficient w_(i). The singularvalue decomposition of Φ is ##EQU9## where σ₁ >σ₂ >. . . >σ_(NL) are thesingular values. The singular value decomposition gives good results,even if the matrix is ill-conditioned. Another advantage of employingthe singular value decomposition is that it allows the use of the methodof principal components.

In systems and control theory literature, the term `principalcomponents` refers to the dominant terms in the singular valuedecomposition of a square matrix, see for example B. C. Moore,"Principal component analysis in linear systems:controllability,observability, and model reduction", IEEE Trans. Autom. Control, Vol.AC-26, pp. 17-31, 1981.

The principal components matrix Φ₁ corresponding to the matrix Φ is##EQU10## where 1≦p≦NL. The principal components matrix Φ₁ provides anoptimal low-rank approximant to Φ. The optimisation problem of eq. [15]can then be solved by determining the inverse Φ₁ ⁻¹.

The use of the principal components approach in combination with thetime delays to define the filter supports allows the filters to focus onwave components (i.e. \ the ground-roll) which have a certain moveoutand coherency over an extended region of the data set. Furthermore,using only a given number of principal components results in only themost dominant coherent components influencing the computation of thefilter coefficients. The number of principal components is an adjustableparameter which allows the user to vary the strength of the filter. Thegreater the number of principal components used, the heavier thefiltering. In general, the smaller the signal-to-noise ratio, theheavier the filtering that should be applied.

An alternative approach of solving the matrix equations is the dampedleast squares solution using: ##EQU11## where I is the NL×NL identitymatrix, and ε² is a small parameter. In the damped least squaresapproach, ill-conditioning is prevented, and the principal componentsresult is approximated, in the sense that components with smallassociated singular values are de-emphasized. The contributions of theprincipal components are only slightly altered. The main advantage ofthis approach is that, since the matrix Φ+ε² I is symmetric, it can besolved using for example the Cholesky decomposition, which is muchfaster than the singular value decomposition.

Whereas in the variant of the invention as described above, themultichannel filter supports are fixed at the beginning of theapplication, depending on the dominant ground-roll velocity, minimumsignal protection gap and filter length. This setup is not ideal for theapplication where the horizontal components are used as reference, sincethe effective moveout between the vertical and horizontal components isnot fixed between the three channels. There it might be advantageous tovary the filter lengths, i.e., to locally adjust the multichannel filtersupports so as to follow the moveout of highest correlation.

In the following, a number of methods for generating the interferencereference channel 30 for use in a method according to the invention willbe described by way of example.

Method 1--Multiple Reference Channels Based on Moveout Differentiation

In a seismic gather, the moveout of the reflected energy in general hasa different moveout than the coherent noise, such as the ground roll.This property can be used to apply adaptive interference cancelling tothe suppression of ground roll. In FIG. 3 an example of traces obtainedfrom seismic survey is shown. The moveout of the seismic reflection 46and the moveout of the ground roll 48 have different gradients. As canbe seen, the ground roll has a much steeper dip than the seismicreflection. On the primary channel 50, i.e. the central trace, theseismic signal is superimposed on a certain phase of the ground roll. Onthe adjacent traces, the same phase of the ground roll does not overlapwith the seismic signal. Thus by appropriate sampling of the adjacenttraces, a noise reference can be obtained for use in adaptiveinterference cancelling. The sampling areas 52, 54 provide support ofthe filter coefficients acting on the reference-channel and are obtainedat a defined distance from the current sample of the primary trace 50.The relative supports of the multichannel filter components can beadjusted by delaying the primary and the reference channel inputs byappropriate amounts. Employing the described method of generating thereference channels, the resulting data-adaptive multichannel filter isshown in FIG. 4, where equivalent elements to those already described inrelation to FIG. 2 are shown with corresponding reference numerals. InFIGS. 3 and 4,

    Δτ.sub.1 =τ.sub.0 -τ.sub.1               [ 16]

and

    Δτ.sub.2 =τ.sub.2 -τ.sub.0 -L            [17]

where L is the number of coefficients used in the adaptive filterassociated with the second sampling.

For illustration purposes, in the rest of this application, the nearesttraces on either side of the trace to be filtered are used as noisereference.

The choice of Δτ_(i) mainly depends on two factors, namely theautocorrelation function of the signature of the seismic reflection, andthe moveout of the ground roll. The autocorrelation function of thesignature of the seismic reflection depends on the power spectrum of theseismic reflection. In general, the window in time defined by Δτ_(i) ischosen to be larger than the lag beyond which the autocorrelationfunction of the signature of the seismic reflection has negligiblevalues. Secondly, the sampling of the nearest traces is chosen so thatthe dominant moveout of the ground roll bisects the sampling area. Thisreduces effects due to the dispersive nature of the ground roll and itsrange of apparent velocities.

Method 2--Multiple Reference Channels Based on Moveout andSpatiotemporal Coherency Differentiation

The target signal and the interference may have both different moveout(apparent velocity) and spatiotemporal coherency. This often occurs insurface seismic acquisition. The ground roll not only has a differentmoveout than the reflection signal, but has significantly morespatiotemporal coherency than the seismic reflections, especially in thetemporal direction. This enables a generalisation to be made fordefining the supports of the multichannel filter components outlined inMethod 1. FIG. 5 shows how the generalisation may be used with tracesobtained from a seismic survey. The moveout of the seismic signal 56 andthe moveout of the ground roll 58 are again of different gradients. Tworegions of filter support are specified for each nearest trace eitherside of the primary channel 60. As for the sampling carried out in FIG.3, the sampling areas, or filter support, 62, 64, 66, 68 leave a band oftime samples around the central sample defined by to. This is so as toprotect the seismic signal.

The ground roll often has temporal coherence over time scales muchgreater than the seismic signal, and in such circumstances appropriateportions 70, 72 of the central trace (ie the primary channel) may beused as a reference.

This generalisation can be very useful in many applications. When tworeference channels are chosen such that one is causal and the otheranti-causal with respect to the incidence direction of the direct groundroll as described in Method 1, the output will contain the reflectingevent, as well as a precursor (a leading event) which is approximatelyhalf the lowpass version of the event, as well as a similar postcursor(a lagging event). If a: minimum-phase response is required at all theprocessing steps, the presence of a precursor is not desirable. If thereference channel is determined by using reference traces which arenearer the source with respect to the primary channel, or by usingsections of reference traces that are causal with respect to the outputpoint, then only a postcursor will be present in the filtered data.

When interference references are generated only by using time delays todefine the filter supports and without the use of any preprocessing, thepresence of reflectors on the reference channels produces smallprecursors and/or postcursors on the processed primary channel.

For example, when two adjacent traces are used for reference a precursorand a postcursor are produced. This effect can also be observed in thefrequency domain: when the reference channels focus on a ground-rollphase with the same polarity as the ground-roll phase in the primarychannel, the low-frequency gain of each filter is approximately -1/K,where K is the number of reference channels. Therefore, there is a nullin the overall response at zero frequency.

The roll-off depends on K; for K=4 the corner frequency is about 5 Hz,for K=2 it's typically about 11-12 Hz. The response also shows some gainin middle frequencies. When the reference channels focus ona ground-rollphase with the opposite polarity as the ground-roll phase in the primarychannel, the low-frequency gain of each filter is approximately +1/K.Then there is no low-frequency roll-off; there is usually a dip inmiddle frequencies.

Since the responses with positive and negative polarity are quitecomplementary, one possibility is to run adaptive interferencecancelling on the same data twice with alternating polarities, and ateach frequency choose the response from the two experiments where themagnitude of the response is the highest, which is greater than unity atmost frequencies. Thus, the deconvolution gain is less than unity atmost frequencies, thus reducing the noise even further while moreideally preserving the seismic reflection signature.

Method 3--Single Reference Channel Based on Median Stack

The median filtering of a sequence of numbers is performed by passing arectangular window over the sequence and replacing each point in thesequence with the median value of the points within the window centeredon the current point. This 1-D median filtering operation can be appliedon a 2-D data set in a specified direction to obtain an estimate of wavecomponents incident in that direction.

The reference channel can be generated by adding the ground roll alongthe dominant moveout using a median stack operation as shown in FIG. 6.The stacking operation is carried out locally on typically 3 to 5neighbouring traces. As shown in FIG. 6, the sensing means, orgeophones, 72, 74, 76 produce reference traces 78 and 80 and a primarytrace 82. The traces are used in a median stack operation. Unlike themean, or conventional stacking operation, the median operation does notsmear the seismic reflections (which do not share the same moveout asthe ground roll), but suppresses them. Thus, the reference channel 90produced by the median stack operation contains a better estimate of theinterfering ground roll with less contamination by the seismic signal ofinterest than is possible using conventional mean stacking. This methodmay be performed on overlapping windows to readjust the stackingmoveout.

Method 4--Multiple Reference Channels Based on Median Stack and MoveoutDifferentiation

The previous embodiment is further developed by firstly applying medianstack 114 and thereafter time shifts τ₁, τ₂ to generate the referencechannels x₁ (t), x₂ (t). FIG. 7 demonstrates how five traces 91, 92, 93,94, 96 can thus be used to generate two reference channels. Signalenergy as leaking through the median filter is further suppressed. Thisembodiment is especially useful when filtering a small number of tracesor data from a small number of geophones or hydrophones.

Method 5--Single Reference Channel Based on Polarisation Differentiation

Ground roll and seismic reflections are often polarised differently,with seismic reflections being linearly polarised and ground roll oftenappearing to be elliptically polarised. Polarisation filtering is usedto distinguish seismic events by the inherent polarisation when theseismic wavefield is measured by multicomponent sensing means, orgeophones.

A method for obtaining a reference channel from polarisation filteringof the components produced from a multicomponent sensing means is shownin FIG. 8. A multicomponent sensing means 116 detects three orthogonalcomponents u_(x) 118, u_(y) 120, u_(z) 122 of the reflected wave. Thepolarisation filter 124 filters the elliptically polarised ground rollsignal and the filtered signal is used as the noise reference 126 inadaptive interference cancelling. The original unfiltered trace 128 isused as the primary channel. Application of the adaptive filter reducesthe residual energy and thus improves ground roll attenuation. Thismethod uses data from a single three component (3 C) geophone and may beapplied to each individual component to recover ground roll attenuationmulticomponent data.

Method 6--Multiple Reference Channels Based on Polarisation and MoveoutDifferentiation

As shown in FIG. 9, this method uses reference geophones 130, 132 of atleast 3 C and a primary geophone 134. The primary geophone produces aprimary channel 136, defined as at time t₀. Polarisation filters 138 and140 are used as in Method 4 to filter the elliptically polarised groundroll signal and produce reference traces 142, 144. The reference traces142, 144 are delayed by t₁ and t₂ so as to sample the reference tracesin accordance with the moveout differentiation method as described inMethod 1. The delayed reference traces are then used to provide thereference channel. It is possible to perform this method using tworeference geophones, the primary channel being obtained from one of thecomponents before filtering.

In this method of obtaining a reference channel, the two completelyindependent criteria of moveout and polarisation are used simultaneouslyto differentiate between the seismic signal and the ground roll. Thismethod'produces a more robust attenuation of ground roll than Method 4.

Method 7--Multiple Reference Channels Generated using Iterated AdaptiveInterference Cancelling

As discussed above, the concept of adaptive interference cancellingrelies on the absence of signal components in the noise referencechannels. In reality, this is rarely the case, there is almost alwayssome signal leakage. The noise estimate from an initial run of adaptiveinterference cancelling may be used as noise reference in a second run,which improves the signal and noise separation obtained from theoriginal data. Use of an iterative approach, using the filtered outputsignal to successively refine the output signal, may be used as manytimes as is deemed advantageous. For methods with multiple referencechannels, this method also has the effect of increasing the effectivearray aperture with only linear increase in computational complexity.

Method 8--Parametric Wavefield Decomposition as Reference Generator

In parametric wavefield decomposition, the problem of decomposing aseismic data set into its constituent wavefields is formulated as aparametric inversion, where each wavefield is modelled by its Fouriercomponents and by frequency-independent parameters (e.g. apparentslowness). The method solves for the wavefield parameters by minimizingthe least squares error between modelled and data waveforms in thefrequency domain. Since the apparent slowness parameters enter theproblem nonlinearly, gradient-based search techniques are used todetermine optimum parameter values. The wavefields are thenreconstructed from the data by solving a linear system at each frequencyfollowed by inverse Fourier transform.

The method of parametric wavefield decomposition can be used toinitially separate the data into its ground-roll and seismic reflectionsignal components. The ground-roll estimate provided by the parametricwavefield decomposition is then used by adaptive interference cancellingas noise reference. Parametric wavefield decomposition in generalachieves good ground-roll and signal separation; however it is affectedby amplitude and time perturbations in the data. Combined use ofparametric wavefield decomposition and adaptive interference cancellingthus gives improved results.

This embodiment can be further developed by firstly applying parametricwavefield decomposition and thereafter time shifts to generate referencechannels. Signal energy that may leak through parametric wavefielddecomposition is further suppressed. The method of parametric wavefielddecomposition is described for example in: C. Esmersoy, "Inversion of Pand SV waves from multicomponent offset vertical seismic profiles,"Geophysics, Vol. 55, pp. 39-50, 1990.

Method 9--Horizontal Components as Reference Channels

In multicomponent recording, most of the compressional body-wave energyis on the vertical component. The horizontal components (transverse andradial) contain ground-roll and shear-wave energy, thus they can beutilized as interference references in adaptive interference cancelling.Unlike the approaches based on polarization filtering which assume thatthe ground-roll is elliptically polarized, this version of 3-C adaptiveinterference cancelling assumes that locally (in the time domain) thereis a certain but unknown relationship between the vertical andhorizontal components of the ground-roll, and effectively estimates thisrelationship. Since seismic reflection signals will also be present inthe horizontal components, the reference channels are generated byapplying time shifts to the horizontal components.

The general idea of using horizontally sensitive geophones assurface-wave detectors to adaptively filter the vertical component isfor example described in a patent by G. A. Crews and D. R. Martinez(U.S. Pat. No. 4,890,264).

The implementation of a method in accordance with the invention will nowbe described for three dimensional surface surveys, i.e. land basedsurveys.

Implementation for 3-D Surveys

In 3-D surface seismic surveys, the current practice is to deploy 2-Daerial arrays of geophones and collect data for varying azimuthaldistributions of source locations. These arrays typically contain of theorder of 24 geophones, and the output signal is obtained by hardwiringthe geophones and summing their outputs.

Adaptive interference cancelling as described in this work may be usedin 3-D surveys, with an aerial pattern of geophones, or net, for examplea hexagon layout of geophones with a central geophone. If the directionof incidence of the direct ground roll (or, equivalently, the sourcelocation) is not supplied at the time of processing, this may beestimated by cross-correlating an arbitrary reference trace with theother traces. Each trace (geophone output) is in turn selected as theprimary channel, and a number of corresponding reference traces arechosen so that sufficient ground roll moveout exists between the primaryand the reference traces. All the filtered traces obtained by using themethod according to the invention are stacked either as in conventionalbeamforming (mean operation), or using median stacking. The delay ofeach filtered trace by different τ₁ and τ₂ values ensures that theprecursor and postcursor features of the signal are in generalde-emphasised, while the seismic reflection is emphasised.

The areal dimensions of these geophones nets are much smaller than thehardwired aerial arrays used at present. In addition the subsequentstacking of the filtered traces preserves more of the signal bandwidththan the current practice, particularly for large source to receiveroffset.

Implementation using Perfect Reconstruction Filter Banks

In most seismic surveys, the coherent noise occupies only a fraction ofthe temporal bandwidth available. For example, in the test data used inthis work, the Nyquist frequency is 250 Hz, while most of the groundroll energy is under 30 Hz. Concentrating filtering efforts to thefrequency band where the coherent noise resides is desirable to reducecomputational cost. One means of achieving this aim involves adding QMF(quadrature mirror filter) perfect reconstruction filter banks to thenoise suppression system using adaptive multichannel filter banks. Thustwo filter banks are effectively involved in the system. The QMF filterbank is used to decompose the traces into frequency bands and decimatebefore adaptive filtering is applied, and then used for resynthesising.Using the QMF filter banks to decimate reduces the number of points tobe processed and also allows reduction in the number of coefficients inthe adaptive filters, bringing in significant savings in computerprocessing time and computer memory requirements.

Due to the frequency overlap needed in the design of QMF filter banks,there is some leakage between bands. The effect of this leakage can beaddressed by widening the band that is processed. Additional suppressionof residual ground roll can be achieved by forming conventional groupswhich can be seen as simple k-domain filters or beamformers. The methodsin accordance with the invention avoid data processing problems wherealiased noise such as ground roll is present.

Due to the frequency overlap needed in the design of QMF filter banks,there is some leakage between bands. The effect of this leakage can beaddressed by widening the band that is processed. Additional suppressionof residual ground roll can be achieved by forming conventional groupswhich can be seen as simple k-domain filters or beamformers. The methodsin accordance with the invention avoid data processing problems wherealiased noise such as ground roll is present.

Aliased ground roll when the geophone (or source) sampling distance isgreater than that required by the Nyquist theorem. In some acquisitiongeometries this may be done by design to reduce the number of geophones(and/or sources) deployed in the field.

Methods which require proper sampling, like k or f-k domain filters arenot able to separate the signal and the ground roll effectively asfrequency domain (k or f-k) filters seek to separate the seismic signaland the ground roll in one step by assigning distinct regions in thetransform domain to the signal and to the noise. The presence ofaliasing prevents this. The method according to the invention does nottransform the data, so it is not affected by aliasing. It focuses on thenoise and removes it; the remaining data is the filtered signal.

The application of the method in accordance with the invention toprocessing of a real data set when compared with processing byconventional methods is now discussed. The data set used was obtainedfrom a common receiver gather, where the spacing between Vibroseis shotpoint locations was 4 meters. In order to be able to inspect the effectof processing on both the signal and the noise before stack, twosynthetic reflections were added. Further evaluation of the processingwas performed after adding arbitrarily spike-type noise to the data set.

Performance of the conventional grouping techniques was simulated byarithmetic averaging of adjacent traces in moving windows. With aconventional group of 6 (20 meter group), some ground roll suppressioncould be achieved, but the distortion of the seismic event hasincreased.

Where the data was analysed using a known normalised LMS algorithm inthe adaptive interference cancellation, the normalised LMS algorithmfailed to track the ground roll. The method of adaptiveinterference-cancelling using the LMS algorithm has limited success insuppressing the ground roll. The LMS algorithm does not converge fastenough to cope with the transient and non-stationary nature of theground roll.

On applying the method according to the present invention to the testdata set, the use of the multichannel filter bank preserves the seismicreflections well, and the ground roll is effectively suppressed. In thisexample K=2, ie there were two reference channels for each primarychannel, and L=30, ie there were 30 coefficients for each channel. Theoutput windows h_(i) (t) were triangular with length 100; Δτ₁ =Δτ₂ =3samples; p=5 was the number of principal components used.

The method in accordance with the invention may be used with incoherentnoise sources other than ground roll, such as steamer noise and rignoise in marine surveys, or sea bottom surveys.

Implementation in Marine Environment

In marine seismic surveys, an acoustic source generates waves whichtravel through the water and into the earth. These are then reflected orrefracted by the sub-surface geological formations, travel back throughthe water and are recorded by long hydrophone arrays which are towednear the surface of the water behind a seismic vessel. The hydrophonesare mounted in streamer cables, or streamers. There are usually 2-12streamers towed which are each several kilometers long. The streamersare made up of 100 meter long sections; each section consists ofhydrophones inside an oil-filled skin. Stress-wires and spacers form theinternal skeleton of the streamer.

While the streamers are being towed behind the vessel, self-noise isgenerated due to a variety of sources. The lurching of the vessel,especially in rough seas, causes vibrations in the stress-wires whichinteract with the connectors and the oil-filled skin, generating bulgewaves (or breathing waves) which propagate down the streamers. Thepressure variations are detected by the hydrophones, adding andcorrupting the detected seismic signals. As the streamer moves throughthe water, boundary layer turbulence causes pressure fluctuations at theouter skin wall, which are again coupled to the hydrophones.

The adaptive interference cancelling methods according to the inventionare directly applicable to the bulge wave. Previous uses of adaptiveinterference cancelling in relation to make seismic surveys have had toutilize special sensors (stress sensors or accelerometers) to form noisereference channels, see U.S. Pat. No. 4,821,241, U.S. Pat. No.5,251,183. The methods according to the invention form noise referencechannels using only the hydrophone data. The methods using 3 componentdata and polarisation filtering are not applicable to the streamer noiseproblem, but the other methods are applicable.

In some marine seismic surveys, there is interference which is incidentin the cross-line direction to the streamers. One such case is thepresence of interference coming from drilling activity on an oil rig.Since there may be as few as two streamers being towed behind thevessel, an array of hydrophones formed in the direction of theinterference would not have sufficient number of hydrophones to utiliseconventional beamforming, or data-independent multichannel filtering,such as f-k filtering. On the other hand, the adaptive interferencecancelling methods according to the invention are directly applicable tothis problem. Along the cross-line array, the interference will havegreater moveout, or apparent slowness than the seismic signals, and thisfeature is exploited in applying the adaptive interference cancellingmethods.

The use of the methods according to the present invention permits alarge reduction in the numbers of geophones or hydrophones required toacquire seismic data.

I claim:
 1. A method for filtering noise from discrete noisy seismicsignals in a time interval ((1, . . . , T)) said method comprising thesteps ofusing said noisy seismic signals for determining at least onereference channel (x(t)) for said time interval as an estimate of saidnoise; determining coefficients for M temporally local filters (w(i,t)),said filters forming a filter bank, and M being a number equal to orlarger than two, by minimizing a cost function (J(t)) representing ameasure of the error of the output of said filter bank; and applyingsaid filter bank to said at least one estimate (x(t)) to determine Mfiltered estimates of said noise.
 2. The method according to claim 1,wherein the cost function (J(t)) is temporally global, representing ameasure of the error of the total output of said filter bank in the timeinterval ((1, . . . , T)).
 3. The method according to claim 1, includingthe step of multiplying M filtered estimates with temporal windowfunctions ((h(i,t)).
 4. The method according to claim 3, wherein thetemporal window functions are characterized by the requirement that onlyadjoining windows overlap.
 5. The method of claim 3, wherein the costfunction (J(t)) is minimized using the approximation that the sum,weighted by window functions, of the output of adjacent filters of the Mfilters (w(i,t)) is equal when applied to the same signal in timeregions where said window functions overlap.
 6. The method of claim 1,wherein the discrete noisy signals used as input are recordings fromindividual seismic sensors prior to any group forming techniques.
 7. Themethod according to claim 1, wherein global optimization of the filteredoutput signal is solved by the method of principal components.
 8. Themethod according to claim 1, wherein the global optimization of thefiltered output signal is solved by a damped least squares approach. 9.The method according to claim 7, wherein the number of principalcomponents is adjusted to vary the degree of filtering.
 10. The methodaccording to claim 1, wherein the M local filters are adaptive.
 11. Themethod according to claim 10, wherein each local filter is amultichannel filter.
 12. The method according to claim 1, wherein thefiltered output signal is reprocessed in an iterative manner to furtherfilter the noise.
 13. The method according to claim 1, wherein thereference signal is obtained by moveout differentiation of data signals.14. The method according to claim 1, wherein the reference signal isobtained by median stacking.
 15. The method according to claim 14,wherein the median stacking is followed by moveout differentiation ofthe stacked signals.
 16. The method according to claim 1, wherein theseismic data is obtained from multicomponent receivers, and thereference signal is obtained by polarisation filtering of each componentsensed by the receivers.
 17. The method according to claim 15, whereinpolarisation filtering is combined with moveout differentiation.